#include <cmath>
#include <limits>
#include <vector>
#include <iostream>

using namespace std;

#include "erf.h"

testerf::erf::erf ()
{
    yeserfc = false;
    largeval = false;
}

long double testerf::erf::geterf(long double z)
{
    long double x, rtn=0.0;
    a = fabs(z);

    x = z * z;
    if (x < 0.875)
    {
        testerf::erf::fill (1);
    }
    else if (x <= 1.0)
    {
        testerf::erf::fill (2);
        a = a - 1.0;
    }
    else if (z < 1.2)
    {
        yeserfc = true;
        testerf::erf::fill (3);
        rtn= testerf::erf::eval (z - 1.0);
        rtn += 1.1609394035130658779364917390740703933002E-5;
        rtn += 0.15728759765625;
        rtn = 1.0 - rtn;
    }
    else if (z < 1.25)
    {
        yeserfc = true;
        testerf::erf::fill (4);
        rtn= testerf::erf::eval (z - 1.125);
        rtn += 8.9850951672359304215530728365232161564636E-6;
        rtn += 0.111602783203125;
        rtn = 1.0 - rtn;
    }
    else
    {
        yeserfc = true;
        largeval = true;

        if((-z * z) < -log(LDBL_MAX))
        {
            cout << "Large value, return value is 2.0 or 0.0.\n";
            if (z < 0)
                return 2.0;
            else
                return 0.0;
        }
        rtn = 1.0 - testerf::erf::largevalue(z);
    }

    if (yeserfc == false)
    {
        int multipler;
        if (z < 0)
            multipler = -1;
        else
            multipler = 1;

        rtn =  multipler * testerf::erf::eval (x);
    }

    testerf::erf::emptyvectors();
    return rtn;
}

long double testerf::erf::eval(long double z)
{
    //TODO: check v1 and v2 for error checks.
    long double y=0.0;
    long double erf_const;
    riter1 = v1.rbegin();
    y = (*riter1);
    ++riter1;
    while (riter1 != v1.rend())
    {
        y = y * z + (*riter1);
        ++riter1;
    }

    long double y2 = 0.0;
    riter1 = v2.rbegin();
    y2 = z* (*riter1);
    ++riter1;
    while (riter1 != v2.rend())
    {
        y2 = y2 * z + (*riter1);
        ++riter1;
    }
    if (largeval)
        return y/y2;
    if (yeserfc)
    {
        return z*y/y2;
    }
    else
    {
        if (z < 0.875)
            erf_const = a;
        else
            erf_const = 0.845062911510467529296875;

        return erf_const + a * y/y2;
    }
}

void testerf::erf::fill(int whattofill)
{
    if (whattofill == 1)
    {
        v1.push_back(-3.858252324254637124543172907442106422373E10);
        v1.push_back(9.580319248590464682316366876952214879858E10);
        v1.push_back(1.302170519734879977595901236693040544854E10);
        v1.push_back(2.922956950426397417800321486727032845006E9);
        v1.push_back(1.764317520783319397868923218385468729799E8);
        v1.push_back(1.573436014601118630105796794840834145120E7);
        v1.push_back(4.028077380105721388745632295157816229289E5);
        v1.push_back(1.644056806467289066852135096352853491530E4);
        v1.push_back(3.390868480059991640235675479463287886081E1);

        v2.push_back(-3.005357030696532927149885530689529032152E11);
        v2.push_back(-1.342602283126282827411658673839982164042E11);
        v2.push_back(-2.777153893355340961288511024443668743399E10);
        v2.push_back(-3.483826391033531996955620074072768276974E9);
        v2.push_back(-2.906321047071299585682722511260895227921E8);
        v2.push_back(-1.653347985722154162439387878512427542691E7);
        v2.push_back(-6.245520581562848778466500301865173123136E5);
        v2.push_back(-1.402124304177498828590239373389110545142E4);
        v2.push_back(-1.209368072473510674493129989468348633579E2);
    }
    else if (whattofill == 2)
    {
        v1.push_back(-4.088889697077485301010486931817357000235E1);
        v1.push_back( 7.157046430681808553842307502826960051036E3);
        v1.push_back(-2.191561912574409865550015485451373731780E3);
        v1.push_back(2.180174916555316874988981177654057337219E3);
        v1.push_back(2.848578658049670668231333682379720943455E2);
        v1.push_back(1.630362490952512836762810462174798925274E2);
        v1.push_back(6.317712353961866974143739396865293596895E0);
        v1.push_back(2.450441034183492434655586496522857578066E1);
        v1.push_back(5.127662277706787664956025545897050896203E-1);

        v2.push_back(1.731026445926834008273768924015161048885E4);
        v2.push_back(1.209682239007990370796112604286048173750E4);
        v2.push_back(1.160950290217993641320602282462976163857E4);
        v2.push_back(5.394294645127126577825507169061355698157E3);
        v2.push_back(2.791239340533632669442158497532521776093E3);
        v2.push_back(8.989365571337319032943005387378993827684E2);
        v2.push_back(2.974016493766349409725385710897298069677E2);
        v2.push_back(6.148192754590376378740261072533527271947E1);
        v2.push_back(1.178502892490738445655468927408440847480E1);
    }
    else if (whattofill == 3)
    {
        v1.push_back(-1.139180936454157193495882956565663294826E3);
        v1.push_back(6.134903129086899737514712477207945973616E2);
        v1.push_back(-4.628909024715329562325555164720732868263E2);
        v1.push_back(4.165702387210732352564932347500364010833E1);
        v1.push_back(-2.286979913515229747204101330405771801610E1);
        v1.push_back(1.870695256449872743066783202326943667722E0);
        v1.push_back(-4.177486601273105752879868187237000032364E0);
        v1.push_back(7.533980372789646140112424811291782526263E-1);
        v1.push_back(-8.629945436917752003058064731308767664446E-2);

        v2.push_back(2.744303447981132701432716278363418643778E3);
        v2.push_back(1.266396359526187065222528050591302171471E3);
        v2.push_back(1.466739461422073351497972255511919814273E3);
        v2.push_back(4.868710570759693955597496520298058147162E2);
        v2.push_back(2.993694301559756046478189634131722579643E2);
        v2.push_back(6.868976819510254139741559102693828237440E1);
        v2.push_back(2.801505816247677193480190483913753613630E1);
        v2.push_back(3.604439909194350263552750347742663954481E0);
    }
    else if (whattofill == 4)
    {
        v1.push_back(-9.652706916457973956366721379612508047640E2);
        v1.push_back(5.577066396050932776683469951773643880634E2);
        v1.push_back(-4.406335508848496713572223098693575485978E2);
        v1.push_back(5.202893466490242733570232680736966655434E1);
        v1.push_back(-1.931311847665757913322495948705563937159E1);
        v1.push_back(-9.364318268748287664267341457164918090611E-2);
        v1.push_back(-3.306390351286352764891355375882586201069E0);
        v1.push_back(7.573806045289044647727613003096916516475E-1);
        v1.push_back(-9.611744011489092894027478899545635991213E-2);

        v2.push_back(3.032829629520142564106649167182428189014E3);
        v2.push_back(1.659648470721967719961167083684972196891E3);
        v2.push_back(1.703545128657284619402511356932569292535E3);
        v2.push_back(6.393465677731598872500200253155257708763E2);
        v2.push_back(3.489131397281030947405287112726059221934E2);
        v2.push_back(8.848641738570783406484348434387611713070E1);
        v2.push_back(3.132269062552392974833215844236160958502E1);
        v2.push_back(4.430131663290563523933419966185230513168E0);
    }
    else if (whattofill == 5)
    {
        v1.push_back(1.293515364043117601705535485785956592493E2);
        v1.push_back(2.474534371269189867053251150063459055230E3);
        v1.push_back(1.756537563959875738809491329250457486510E4);
        v1.push_back(5.977479535376639763773153344676726091607E4);
        v1.push_back(1.054313816139671870123172936972055385049E5);
        v1.push_back(9.754699773487726957401038094714603033904E4);
        v1.push_back(4.579131242577671038339922925213209214880E4);
        v1.push_back(1.000710322164510887997115157797717324370E4);
        v1.push_back(8.496863250712471449526805271633794700452E2);
        v1.push_back(1.797349831386892396933210199236530557333E1);

        v2.push_back(2.292696320307033494820399866075534515002E2);
        v2.push_back(4.500632608295626968062258401895610053116E3);
        v2.push_back(3.321218723485498111535866988511716659339E4);
        v2.push_back(1.196084512221845156596781258490840961462E5);
        v2.push_back(2.287033883912529843927983406878910939930E5);
        v2.push_back(2.370223495794642027268482075021298394425E5);
        v2.push_back(1.305173734022437154610938308944995159199E5);
        v2.push_back(3.589386258485887630236490009835928559621E4);
        v2.push_back(4.339996864041074149726360516336773136101E3);
        v2.push_back(1.753135522665469574605384979152863899099E2);
    }
}

testerf::erf::~erf()
{
    testerf::erf::emptyvectors();
}

inline void testerf::erf::emptyvectors()
{
    if (v1.empty() == false)
        v1.clear();
    if (v2.empty() == false)
        v2.clear();
}

long double testerf::erf::largevalue(long double x)
{
    long double p,y,z;
    z = 1.0/(x*x);
    testerf::erf::fill (5);
    p = testerf::erf::eval (z);
    y = p / x;
    long double u, u1, m, f;
    const long double M = 32768.0;
    const long double MINV = 3.0517578125e-5;
    x = fabsl (x);
    x = -x;
    m = MINV * floorl(M * x + 0.5);
    f = x - m;
    u = m * m;
    u1 = 2 * m * f  +  f * f;
    u = -u;
    u1 = -u1;
    u = exp(u) * exp(u1);
    return(y * u);
}
